A review of statistical computing with TMB

# A review of statistical computing with TMB
### Andrea Havron<br>NOAA Fisheries, OST

---

.footnote[U.S. Department of Commerce | National Oceanic and Atmospheric Administration | National Marine Fisheries Service]

code.cpp{
  font-size: 14px;
}
code.r{
  font-size: 14px;
}

</style>

---
# Statistical Computing Review
<br>
.pull-left[
<img src="data:image/png;base64,#static/TMB-trifecta.png" width="60%" style="display: block; margin: auto auto auto 0;" />
]

---

### Review:

---
# ML Inference

#### What is the likelihood of getting 30 heads in 100 coin flips?

#### What is the likelihood of getting 30 heads in 100 coin flips?

.pull-left[
1. Specify the model
2. **Calculate the likelihood**<br><br>
`$L(p; n, y) = \frac{n!}{y!(n-y)!}p^y(1-p)^{n-y}$`
]

.pull-right[
`$y = dbinom(x = 30, n = 100, p)$`
<div id="htmlwidget-1c498dc220e57b77b6c0" style="width:450px;height:375px;" class="plotly html-widget"></div>
<script type="application/json" data-for="htmlwidget-1c498dc220e57b77b6c0">{"x":{"visdat":{"48b472a20d2":["function () ","plotlyVisDat"]},"cur_data":"48b472a20d2","attrs":{"48b472a20d2":{"x":{},"y":{},"mode":"lines","alpha_stroke":1,"sizes":[10,100],"spans":[1,20],"type":"scatter"}},"layout":{"width":450,"height":375,"margin":{"b":40,"l":60,"t":25,"r":10},"xaxis":{"domain":[0,1],"automargin":true,"title":"p"},"yaxis":{"domain":[0,1],"automargin":true,"title":"y"},"hovermode":"closest","showlegend":false},"source":"A","config":{"modeBarButtonsToAdd":["hoverclosest","hovercompare"],"showSendToCloud":false},"data":[{"x":[0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,0.2,0.21,0.22,0.23,0.24,0.25,0.26,0.27,0.28,0.29,0.3,0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39,0.4,0.41,0.42,0.43,0.44,0.45,0.46,0.47,0.48,0.49,0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.61,0.62,0.63,0.64,0.65,0.66,0.67,0.68,0.69,0.7,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79,0.8,0.81,0.82,0.83,0.84,0.85,0.86,0.87,0.88,0.89,0.9,0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99,1],"y":[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.001,0.001,0.003,0.005,0.009,0.015,0.024,0.034,0.046,0.058,0.069,0.079,0.085,0.087,0.085,0.079,0.071,0.06,0.049,0.039,0.029,0.021,0.015,0.01,0.006,0.004,0.002,0.001,0.001,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],"mode":"lines","type":"scatter","marker":{"color":"rgba(31,119,180,1)","line":{"color":"rgba(31,119,180,1)"}},"error_y":{"color":"rgba(31,119,180,1)"},"error_x":{"color":"rgba(31,119,180,1)"},"line":{"color":"rgba(31,119,180,1)"},"xaxis":"x","yaxis":"y","frame":null}],"highlight":{"on":"plotly_click","persistent":false,"dynamic":false,"selectize":false,"opacityDim":0.2,"selected":{"opacity":1},"debounce":0},"shinyEvents":["plotly_hover","plotly_click","plotly_selected","plotly_relayout","plotly_brushed","plotly_brushing","plotly_clickannotation","plotly_doubleclick","plotly_deselect","plotly_afterplot","plotly_sunburstclick"],"base_url":"https://plot.ly"},"evals":[],"jsHooks":[]}</script>
]
---

# ML Inference

**What is the likelihood of getting 30 heads in 100 coin flips?**

.pull-left[
1. Specify the model
2. Calculate the likelihood
3. **Calculate the negative log-likelihood**<br><br>
`$-\ell(p; n, y) = -[ln\big(\frac{n!}{y!(n-y)!}\big) + yln(p)$`<br>
     `$$+ (n-y)ln(1-p)]$$`
]

.pull-right[
`$nll = -log(dbinom(x = 30, n = 100, p))$`
<div id="htmlwidget-eb987958ea719f25d5b4" style="width:450px;height:375px;" class="plotly html-widget"></div>
<script type="application/json" data-for="htmlwidget-eb987958ea719f25d5b4">{"x":{"visdat":{"48b47dde6521":["function () ","plotlyVisDat"]},"cur_data":"48b47dde6521","attrs":{"48b47dde6521":{"x":{},"y":{},"mode":"lines","alpha_stroke":1,"sizes":[10,100],"spans":[1,20],"type":"scatter"}},"layout":{"width":450,"height":375,"margin":{"b":40,"l":60,"t":25,"r":10},"xaxis":{"domain":[0,1],"automargin":true,"title":"p"},"yaxis":{"domain":[0,1],"automargin":true,"title":"nll"},"hovermode":"closest","showlegend":false},"source":"A","config":{"modeBarButtonsToAdd":["hoverclosest","hovercompare"],"showSendToCloud":false},"data":[{"x":[0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,0.2,0.21,0.22,0.23,0.24,0.25,0.26,0.27,0.28,0.29,0.3,0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39,0.4,0.41,0.42,0.43,0.44,0.45,0.46,0.47,0.48,0.49,0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.61,0.62,0.63,0.64,0.65,0.66,0.67,0.68,0.69,0.7,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79,0.8,0.81,0.82,0.83,0.84,0.85,0.86,0.87,0.88,0.89,0.9,0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99,1],"y":[null,80.217,60.133,48.687,40.782,34.82,30.092,26.216,22.966,20.198,17.811,15.734,13.914,12.313,10.899,9.648,8.54,7.56,6.693,5.93,5.261,4.678,4.174,3.744,3.382,3.084,2.847,2.668,2.542,2.468,2.444,2.468,2.537,2.651,2.808,3.007,3.248,3.528,3.848,4.207,4.604,5.04,5.514,6.025,6.575,7.162,7.787,8.45,9.152,9.893,10.673,11.493,12.354,13.256,14.201,15.189,16.221,17.299,18.425,19.599,20.823,22.099,23.43,24.817,26.262,27.769,29.34,30.979,32.688,34.473,36.336,38.284,40.321,42.453,44.686,47.029,49.489,52.076,54.801,57.675,60.713,63.931,67.347,70.985,74.869,79.032,83.51,88.351,93.611,99.363,105.7,112.743,120.66,129.683,140.153,152.598,167.904,187.731,215.806,264.021,null],"mode":"lines","type":"scatter","marker":{"color":"rgba(31,119,180,1)","line":{"color":"rgba(31,119,180,1)"}},"error_y":{"color":"rgba(31,119,180,1)"},"error_x":{"color":"rgba(31,119,180,1)"},"line":{"color":"rgba(31,119,180,1)"},"xaxis":"x","yaxis":"y","frame":null}],"highlight":{"on":"plotly_click","persistent":false,"dynamic":false,"selectize":false,"opacityDim":0.2,"selected":{"opacity":1},"debounce":0},"shinyEvents":["plotly_hover","plotly_click","plotly_selected","plotly_relayout","plotly_brushed","plotly_brushing","plotly_clickannotation","plotly_doubleclick","plotly_deselect","plotly_afterplot","plotly_sunburstclick"],"base_url":"https://plot.ly"},"evals":[],"jsHooks":[]}</script>
]
---

# ML Inference

What is the likelihood of getting 30 heads in 100 coin flips?

.pull-left[
1. Specify the model
2. Calculate the likelihood
3. Calculate the negative log-likelihood
`$-\ell(p; n, y) = -[ln\big(\frac{n!}{y!(n-y)!}\big) + yln(p)$`<br>
     `$$+ (n-y)ln(1-p)]$$`
4. **Calculate the derivative w.r.t. `$p$`**<br><br>
`$\frac{d(\ell(p; n, y))}{dp} = \frac{y}{p}- \frac{n-y}{1-p}$`
]

.pull-right[
`$-ln\big[L(p; n = 100, y = 30)\big]$`
<img src="data:image/png;base64,#00_00_Statistical_Computing_Review_files/figure-html/unnamed-chunk-4-1.png" width="80%" />
]
---

# ML Inference

What is the likelihood of getting 30 heads in 100 coin flips?

.pull-left[
1. Specify the model
2. Calculate the likelihood
3. Calculate the negative log-likelihood
4. Calculate the derivate wrt p
5. **Set to 0 and solve for MLE**<br>
`$0 = \frac{y}{p}- \frac{n-y}{1-p}$` <br>
`$E[p] = \frac{y}{n}$`<br>
`$E[y] = np$`
]

.pull-right[
`$-ln\big[L(p; N = 100, y = 30)\big]$`
<img src="data:image/png;base64,#00_00_Statistical_Computing_Review_files/figure-html/unnamed-chunk-5-1.png" width="70%" />
<br>
`$\hat{p} = \frac{30}{100} = 0.3$`
]

---

# ML Inference

What is the likelihood of getting 30 heads in 100 coin flips?

.pull-left[
1. Specify the model
2. Calculate the likelihood
3. Calculate the negative log-likelihood
4. Calculate the derivate wrt p
5. Set to 0 and solve for MLE<br>
`$0 = \frac{y}{p}- \frac{n-y}{1-p}$` <br>
`$E[y] = np$`
6. **Approximate Var[p] using the second derivative**<br>
`$-\frac{y}{p^2} - \frac{(n-y)}{(1-p)^2}$`<br>
`$-\frac{np}{p^2} - \frac{(n-np)}{(1-p)^2}$`<br>
`$-\frac{n}{p} - \frac{n(1-p)}{1-p}$`<br>
`$l''(p) = -\frac{n(1-p)}{p}$`<br>
`$Var[p] = \frac{p}{n(1-p)}$`
]

.pull-right[
`$Var[p] \approx -\frac{1}{l''(p)}$`<br>
`$SE[p] \approx \sqrt{ \frac{.3}{100(.7)}} \approx 0.065$`<br>
<img src="data:image/png;base64,#00_00_Statistical_Computing_Review_files/figure-html/unnamed-chunk-6-1.png" width="70%" />
]
---
#Multivariate asymptotics

* For N-d models, the curvature is represented by a NxN **Hessian** matrix of 2nd partial derivatives
* Inverting the negative Hessian gives us a covariance matrix

.pull-left[
`\begin{align}
(\mathbb{H}_{f})_{i,j} &= \frac{\partial^2f}{\partial \theta_{i}, \partial x\theta_{j}} = \frac{-1}{Var(\Theta)}
\end{align}`
<br>
[**What causes a singular covariance matrix?**](https://andrea-havron.shinyapps.io/mvnplot/)
]

.pull-right[
![](data:image/png;base64,#00_00_Statistical_Computing_Review_files/figure-html/mvnorm-1.png)
]

---
# What is AD?

.three-column-left[
**Automatic Differentiation**<br>
Derivatives calculated automatically using the chain rule<br>
.p[
- Efficient: forward mode, O(n); reverse mode, O(m)
- Accurate 
- Higher order derivatives: easy
]]
.three-column-left[
**Symbolic Differentiation**<br>
Computer program converts function into exact derivative function<br>
.p[
- Inefficient: O(2n) for n input variables
- Exact 
- Higher order derivatives: difficult due to complexity
]]
.three-column-left[
**Numerical Differentiation**<br>
Approximation that relies on finite differences
.p[
- Efficient: O(n) 
- Trade-off between truncation error versus round-off error
- Higher order derivatives calculation difficult due to error accumulation
]]

---
# Computational Graph (Tape)
<br>
.three-column-left[

```cpp
//Program
v1: x = ?
v2: y = ?
v3: a = x * y
v4: b = sin(y)
v5: z = a + b
```
]
.three-column[
<img src="data:image/png;base64,#static/comp-graph.png" width="100%" style="display: block; margin: auto auto auto 0;" />
]
.three-column-left[

```cpp
//Reverse Mode
dz = ?
da = dz
db = dz
dx = yda
dy = xda + cos(y)db
```
]

---
# Reverse Mode
.pull-left[
**Static (TMB: CppAD, TMBad)**<br>
The graph is constructed once before execution
.p[
- Less flexibility with conditional statements that depend on parameters. 
- Atomic functions can be used when conditional statements depend on parameters
- High portability 
- Graph optimization possible
]]

.pull-right[
**Dynamic (Stan: Stan Math Library, ADMB: AUTODIF)**<br>
The graph is defined as forward computation is executed at every iteration
.p[
- Flexibility with conditional statements
- Optimization routine implemented into executable
- Less time for graph optimization
]
]

---

# Type Systems in R and C++
---
# Dynamic vs. Static Typing
<br>
.pull-left[
**R: Dynamic**
.p[
- Type checking occurs at run time
- The values and types associated with names can change
- Change in type tends to be implicit
]]

.pull-right[
**C++: Static**
.p[
- Type checking occurs at compile time
- The values associated with a given name can be limited to just a few types and may be unchangeable
- Change in type tends to be explicit
]]

---
# What is Templated C++?
<br><br>
* Generic programming
* Allows developer to write functions and classes that are independent of Type
* Templates are expanded at compile time

```cpp

template <class T>
  T add(T x, T y){
  return x + y;
}

int main(){
  int a = 1;
  int b = 2;
  double c = 1.1;
  double d = 2.1;
  int d = add(a,b);
  double e = add(c,d);
}

```
]
.pull-right[

```cpp
int add(int x, int y){
  return x + y;
}
double add(double x, double y){
  return x + y;
}
```
]

---

# TMB AD Systems
<br>
.pull-left[
**CppAD**
- [CppAD package](https://coin-or.github.io/CppAD/doc/cppad.htm)
]
.pull-right[
**TMBad**<br>
- TMBad is available with TMB 1.8.0 and higher
]